Computer Graphics Methods and Systems Using Quasi-Monte Carlo Methodology

ABSTRACT

A computer graphics system generates a pixel value for a pixel in an image, the pixel value being representative of a point in a scene as recorded on an image plane of a simulated camera, the computer graphics system comprising a sample point generator and a function evaluator. The sample point generator is configured to generate a set of sample points, at least one sample point being generated using at least one dependent sample comprising at least one element of a low-discrepancy sequence offset by at least one element of another low-discrepancy sequence. The function evaluator is configured to generate at least one value representing an evaluation of a selection function at one of the sample points, the value generated by the function evaluator corresponding to the pixel value.

CROSS-REFERENCE TO RELATED APPLICATIONS

This application is a Continuation of U.S. patent application Ser. No.11/619,275 filed Jan. 3, 2007, which is a Continuation of co-pendingU.S. patent application Ser. No. 10/299,958 filed Nov. 19, 2002, whichis a Continuation-in-Part of U.S. patent application Ser. No. 09/884,861filed Jun. 19, 2001, which claims priority from U.S. Provisional PatentApps. Ser. Nos. 60/265,934 filed Feb. 1, 2001 and 60/212,286 filed Jun.19, 2000, each and all of which are incorporated by reference herein asif set forth in their entireties. Also incorporated by reference hereinas if set forth in its entirety is U.S. patent application Ser. No.08/880,418 filed Jun. 23, 1997 (hereinafter “the Grabensteinapplication”), assigned to the assignee of this application.

FIELD OF THE INVENTION

The invention relates generally to computer graphics, and moreparticularly, to computer graphics systems, methods, devices andcomputer program products that generate pixel values, for generatingimages, using quasi-Monte Carlo methodologies. One aspect of theinvention provides a methodology that makes use of trajectory splittingby dependent sampling.

BACKGROUND OF THE INVENTION

In computer graphics, a computer is used to generate digital data thatrepresents the projection of surfaces of objects in, for example, athree-dimensional scene, illuminated by one or more light sources, ontoa two-dimensional image plane, to simulate the recording of the sceneby, for example, a camera. The camera may include a lens for projectingthe image of the scene onto the image plane, or it may comprise apinhole camera in which case no lens is used. The two-dimensional imageis in the form of an array of picture elements (“pixels” or “pels”), andthe digital data generated for each pixel represents the color andluminance of the scene as projected onto the image plane at the point ofthe respective pixel in the image plane. The surfaces of the objects mayhave any of a number of characteristics, including shape, color,specularity, texture, and so forth, which are preferably rendered in theimage as closely as possible, to provide a realistic-looking image.

Generally, the contributions of the light reflected from the variouspoints in the scene to the pixel value representing the color andintensity of a particular pixel are expressed in the form of the one ormore integrals of relatively complicated functions. Since the integralsused in computer graphics generally will not have a closed-formsolution, numerical methods must be used to evaluate them and therebygenerate the pixel value. Typically, a conventional “Monte Carlo” methodhas been used in computer graphics to numerically evaluate theintegrals. Generally, in the Monte Carlo method, to evaluate an integral

$\begin{matrix}{{\langle f\rangle} = {\int_{{({0,1})}^{s}}{{f(x)}{x}}}} & (1)\end{matrix}$

where f(x) is a real function on the “s”-dimensional unit cube [0,1)^(s)(that is, an s-dimensional cube each of whose dimension includes “zero,”and excludes “one”), first a number “N” statistically-independentrandomly-positioned points xi, i=1, . . . , N, are generated over theintegration domain. The random points xi are used as sample points forwhich sample values f(xi) are generated for the function f(x), and anestimate ƒ for the integral is generated as

$\begin{matrix}{{{{\langle f\rangle} \approx \hat{f}} = {\frac{1}{N}{\overset{N}{\sum\limits_{i = 1}}{f\left( x_{i} \right)}}}}} & (2)\end{matrix}$

As the number of random points used in generating the sample pointsf(xi) increases, the value of the estimate ƒ will converge toward theactual value of the integral

ƒ

. Generally, the distribution of estimate values that will be generatedfor various values of “N,” that is, for various numbers of samplepoints, of being normal distributed around the actual value with astandard deviation σ which can be estimated by

$\begin{matrix}{\sigma = \sqrt{\frac{1}{N - 1}\left( {\overset{\_}{f^{2}} - {\overset{\_}{f}}^{2}} \right)}} & (3)\end{matrix}$

if the points x_(i) used to generate the sample values f(x_(i)) arestatistically independent, that is, if the points x_(i) are trulypositioned at random in the integration domain.

Generally, it has been believed that random methodologies like the MonteCarlo method are necessary to ensure that undesirable artifacts, such asMoiré patterns and aliasing and the like, which are not in the scene,will not be generated in the generated image. However, several problemsarise from use of the Monte Carlo method in computer graphics. First,since the sample points x; used in the Monte Carlo method are randomlydistributed, they may clump in various regions over the domain overwhich the integral is to be evaluated. Accordingly, depending on the setof points that are generated, in the Monte Carlo method for significantportions of the domain there may be no sample points x_(i) for whichsample values f(x_(i)) are generated. In that case, the error can becomequite large. In the context of generating a pixel value in computergraphics, the pixel value that is actually generated using the MonteCarlo method may not reflect some elements which might otherwise bereflected if the sample-points x_(i) were guaranteed to be more evenlydistributed over the domain. This problem can be alleviated somewhat bydividing the domain into a plurality of sub-domains, but it is generallydifficult to determine a priori the number of sub-domains into which thedomain should be divided, and, in addition, in a multi-dimensionalintegration region, which would actually be used in computer graphicsrendering operations, the partitioning of the integration domain intosub-domains, which are preferably of equal size, can be quitecomplicated.

In addition, since the method makes use of random numbers, the error |ƒ−

ƒ

| (where |x| represents the absolute value of the value “x”) between theestimate value ƒ and actual value

ƒ

f) is probabilistic, and, since the error values for various largevalues of “N” are close to normal distribution around the actual value

ƒ

, only sixty-eight percent of the estimate values ƒ that might begenerated are guaranteed to lie within one standard deviation of theactual value

ƒ

.

Furthermore, as is clear from equation (3), the standard deviation σdecreases with increasing numbers “N” of sample points, proportional tothe reciprocal of square root of “N” (that is, 1/√{square root over(N)}). Thus, if it is desired to reduce the statistical error by afactor of two, it will be necessary to increase the number of samplepoints N by a factor of four, which, in turn, increases thecomputational load that is required to generate the pixel values, foreach of the numerous pixels in the image.

Additionally, since the Monte Carlo method requires random numbers todefine the coordinates of respective sample points x_(i) in theintegration domain, an efficient mechanism for generating random numbersis needed. Generally, digital computers are provided with so-called“random number” generators, which are computer programs which can beprocessed to generate a set of numbers that are approximately random.Since the random number generators use deterministic techniques, thenumbers that are generated are not truly random. However, the propertythat subsequent random numbers from a random number generator arestatistically independent should be maintained by deterministicimplementations of pseudo-random numbers on a computer.

The Grabenstein application describes a computer graphics system andmethod for generating pixel values for pixels in an image using astrictly deterministic methodology for generating sample points, whichavoids the above-described problems with the Monte Carlo method. Thestrictly deterministic methodology described in the Grabensteinapplication provides a low-discrepancy sample point sequence whichensures, a priori, that the sample points are generally more evenlydistributed throughout the region over which the respective integralsare being evaluated. In one embodiment, the sample points that are usedare based on a so-called Halton sequence. See, for example, J. H.Halton, Numerische Mathematik, Vol. 2, pp. 84-90 (1960) and W. H. Press,et al., Numerical Recipes in Fortran (2d Edition) page 300 (CambridgeUniversity Press, 1992). In a Halton sequence generated for number base“b,” where base “b” is a selected prime number, the “k-th” value of thesequence, represented by H_(b) ^(k) is generated by use of a “radicalinverse” function radical inverse function Φ_(b) that is generallydefined as

$\begin{matrix}{{{\Phi_{b}\text{:}N_{0}}->I}i = \left. {\sum\limits_{j = 0}^{\infty}{{a_{j}(i)}b^{j}}}\mapsto{\sum\limits_{j = 0}^{\infty}{{a_{j}(i)}b^{{- j} - 1}}} \right.} & (4)\end{matrix}$

where (α_(j))_(j=0) ^(∞) is the representation of “i” in integer base“b.” Generally, a radical inverse of a value “k” is generated by

-   -   (1) writing the value “k” as a numerical representation of the        value in the selected base “b,” thereby to provide a        representation for the value as D_(M)D_(M−1) . . . D₂D₁, where        D_(m) (m=1, 2, . . . , M) are the digits of the representation,    -   (2) putting a radix point (corresponding to a decimal point for        numbers written in base ten) at the least significant end of the        representation D_(M)D_(M−1) . . . D₂D₁ written in step (1)        above, and    -   (3) reflecting the digits around the radix point to provide        0.D₁D₂ . . . D_(M−1)D_(M), which corresponds to H_(b) ^(k).        It will be appreciated that, regardless of the base “b” selected        for the representation, for any series of values, one, two . . .        “k,” written in base “b,” the least significant digits of the        representation will change at a faster rate than the most        significant digits. As a result, in the Halton sequence H_(b) ¹,        H_(b) ², . . . H_(b) ^(k), the most significant digits will        change at the faster rate, so that the early values in the        sequence will be generally widely distributed over the interval        from zero to one, and later values in the sequence will fill in        interstices among the earlier values in the sequence. Unlike the        random or pseudo-random numbers used in the Monte Carlo method        as described above, the values of the Halton sequence are not        statistically independent; on the contrary, the values of the        Halton sequence are strictly deterministic, “maximally avoiding”        each other over the interval, and so they will not clump,        whereas the random or pseudo-random numbers used in the Monte        Carlo method may clump.

It will be appreciated that the Halton sequence as described aboveprovides a sequence of values over the interval from zero to one,inclusive along a single dimension. A multi dimensional Halton sequencecan be generated in a similar manner, but using a different base foreach dimension, where the bases are relatively prime.

A generalized Halton sequence, of which the Halton sequence describedabove is a special case, is generated as follows. For each startingpoint along the numerical interval from zero to one, inclusive, adifferent Halton sequence is generated. Defining the pseudo-sum x⊕_(p)yfor any x and y over the interval from zero to one, inclusive, for anyinteger “p” having a value greater than two, the pseudo-sum is formed byadding the digits representing “x” and “y” in reverse order, from themost-significant digit to the least-significant digit, and for eachaddition also adding in the carry generated from the sum of next moresignificant digits. Thus, if “x” in base “b” is represented by 0.X₁X₂ .. . X_(M−1)X_(M), where each “X_(m)” is a digit in base “b,” and if “y”in base “b” is represented by 0.Y₁Y₂ . . . Y_(N−1)Y_(N), where each“Y_(n)” is a digit in base “b” (and where “M,” the number of digits inthe representation of “x” in base “b”, and “N,” the number of digits inthe representation of “y” in base “b”, may differ), then the pseudo-sum“z” is represented by 0.Z₁Z₂ . . . Z_(L−1)Z_(L), where each “Z₁” is adigit in base “b” given by Z₁=(X₁+Y₁+C₁) mod b, where “mod” representsthe modulo function, and

$C_{1} = \left\{ \begin{matrix}1 & {{{{for}\mspace{14mu} X_{l - 1}} + Y_{l - 1} + Z_{l - 1}} \geq b} \\0 & {otherwise}\end{matrix} \right.$

is a carry value from the “l−1st” digit position, with C₁ being set tozero.

Using the pseudo-sum function as described above, the generalized Haltonsequence that is used in the system described in the Grabensteinapplication is generated as follows. If “b” is an integer, and x₀ is anarbitrary value on the interval from zero to one, inclusive, then the“p”-adic von Neumann-Kakutani transformation T_(b)(x) is given by

$\begin{matrix}{{{T_{p}(x)}:={x \oplus_{p}\frac{1}{b}}},} & (5)\end{matrix}$

and the generalized Halton sequence x₀, x₁, x₂, . . . is definedrecursively as

x _(n+1) =T _(b)(x _(n))  (6)

From equations (5) and (6), it is clear that, for any value for “b,” thegeneralized Halton sequence can provide that a different sequence willbe generated for each starting value of “x,” that is, for each x₀. Itwill be appreciated that the Halton sequence H_(b) ^(k) as describedabove is a special case of the generalized Halton sequence (equations(5) and (6)) for x₀=0.

The use of a strictly deterministic low-discrepancy sequence such as theHalton sequence or the generalized Halton sequence can provide a numberof advantages over the random or pseudo-random numbers that are used inconnection with the Monte Carlo technique. Unlike the random numbersused in connection with the Monte Carlo technique, the low discrepancysequences ensure that the sample points are more evenly distributed overa respective region or time interval thereby reducing error in the imagewhich can result from clumping of such sample points which can occur inthe Monte Carlo technique. That can facilitate the generation of imagesof improved quality when using the same number of sample points at thesame computational cost as in the Monte Carlo technique.

SUMMARY OF THE INVENTION

The invention provides a new and improved system andcomputer-implemented method for evaluating integrals using a quasi-MonteCarlo methodology that makes use of trajectory splitting by dependentsampling.

In brief summary, the invention provides a computer graphics system forgenerating a pixel value for a pixel in an image, the pixel value beingrepresentative of a point in a scene as recorded on an image plane of asimulated camera, the computer graphics system comprising a sample pointgenerator and a function evaluator. The sample point generator isconfigured to generate a set of sample points, at least one sample pointbeing generated using at least one dependent sample, the at least onedependent sample comprising at least one element of a low-discrepancysequence offset by at least one element of another low-discrepancysequence. The function evaluator is configured to generate at least onevalue representing an evaluation of a selected function at one of thesample points generated by the sample point generator, the valuegenerated by the function evaluator corresponding to the pixel value.

BRIEF DESCRIPTION OF THE DRAWINGS

This invention is pointed out with particularity in the appended claims.The above and further advantages of this invention may be betterunderstood by referring to the following description taken inconjunction with the accompanying drawings, in which:

FIG. 1 depicts an illustrative computer graphics system that evaluatesintegrals using a quasi-Monte Carlo methodology that makes use oftrajectory splitting by dependent sampling.

DETAILED DESCRIPTION OF AN ILLUSTRATIVE EMBODIMENT

The invention provides an computer graphic system and method forgenerating pixel values for pixels in an image of a scene, which makesuse of a strictly-deterministic quasi-Monte Carlo methodology that makesuse of trajectory splitting by dependent sampling for generating samplepoints for use in generating sample values for evaluating the integralor integrals whose function(s) represent the contributions of the lightreflected from the various points in the scene to the respective pixelvalue, rather than the random or pseudo-random Monte Carlo methodologywhich has been used in the past. The strictly-deterministic methodologyensures apiori that the sample points will be generally more evenlydistributed over the interval or region over which the integral(s) is(are) to be evaluated in a low-discrepancy manner.

FIG. 1 attached hereto depicts an illustrative computer system 10 thatmakes use of such a strictly deterministic methodology. With referenceto FIG. 1, the computer system 10 in one embodiment includes a processormodule 11 and operator interface elements comprising operator inputcomponents such as a keyboard 12A and/or a mouse 12B (generallyidentified as operator input element(s) 12) and an operator outputelement such as a video display device 13. The illustrative computersystem 10 is of the conventional stored-program computer architecture.The processor module 11 includes, for example, one or more processor,memory and mass storage devices, such as disk and/or tape storageelements (not separately shown), which perform processing and storageoperations in connection with digital data provided thereto. Theoperator input element(s) 12 are provided to permit an operator to inputinformation for processing. The video display device 13 is provided todisplay output information generated by the processor module 11 on ascreen 14 to the operator, including data that the operator may inputfor processing, information that the operator may input to controlprocessing, as well as information generated during processing. Theprocessor module 11 generates information for display by the videodisplay device 13 using a so-called “graphical user interface” (“GUI”),in which information for various applications programs is displayedusing various “windows.” Although the computer system 10 is shown ascomprising particular components, such as the keyboard 12A and mouse 12Bfor receiving input information from an operator, and a video displaydevice 13 for displaying output information to the operator, it will beappreciated that the computer system 10 may include a variety ofcomponents in addition to or instead of those depicted in FIG. 1.

In addition, the processor module 11 includes one or more network ports,generally identified by reference numeral 14, which are connected tocommunication links which connect the computer system 10 in a computernetwork. The network ports enable the computer system 10 to transmitinformation to, and receive information from, other computer systems andother devices in the network. In a typical network organized accordingto, for example, the client-server paradigm, certain computer systems inthe network are designated as servers, which store data and programs(generally, “information”) for processing by the other, client computersystems, thereby to enable the client computer systems to convenientlyshare the information. A client computer system which needs access toinformation maintained by a particular server will enable the server todownload the information to it over the network. After processing thedata, the client computer system may also return the processed data tothe server for storage. In addition to computer systems (including theabove-described servers and clients), a network may also include, forexample, printers and facsimile devices, digital audio or video storageand distribution devices, and the like, which may be shared among thevarious computer systems connected in the network. The communicationlinks interconnecting the computer systems in the network may, as isconventional, comprise any convenient information-carrying medium,including wires, optical fibers or other media for carrying signalsamong the computer systems. Computer systems transfer information overthe network by means of messages transferred over the communicationlinks, with each message including information and an identifieridentifying the device to receive the message.

It will be helpful to initially provide some background on operationsperformed by the computer graphics system in generating an image.Generally, the computer graphic system generates an image that attemptsto simulate an image of a scene that would be generated by a camera. Thecamera includes a shutter that will be open for a predetermined time Tstarting at a time t₀ to allow light from the scene to be directed to animage plane. The camera may also include a lens or lens model(generally, “lens”) that serves to focus light from the scene onto theimage plane. The average radiance flux L_(m,n) through a pixel atposition (m,n) on an image plane P, which represents the plane of thecamera's recording medium, is determined by

$\begin{matrix}{L_{m,n} = {\frac{1}{{A_{P}} \cdot T \cdot {A_{L}}}{\int_{A_{P}}{\overset{t_{0} + T}{\int\limits_{t_{0}}}{\int_{A_{L}}{{L\left( {{h\left( {x,t,y} \right)},{- {\omega \left( {x,t,y} \right)}}} \right)}{f_{m,n}\left( {x,y,t} \right)}{y}{t}{x}}}}}}} & (7)\end{matrix}$

where “A_(P)” refers to the area of the pixel, A_(L) refers to the areaof the portion of the lens through which rays of light pass from thescene to the pixel, and f_(m,n) represents a filtering kernel associatedwith the pixel. An examination of the integral in equation (7) willreveal that, for the variables of integration, “x,” “y” and “t,” thevariable “y” refers to integration over the lens area (A_(L)), thevariable “t” refers to integration over time (the time interval from tot₀ t₀+T) and the variable “x” refers to integration over the pixel area(A_(P)).

The value of the integral in equation (7) is approximated in accordancewith a quasi-Monte Carlo methodology by identifying N_(P) sample pointsx_(i) in the pixel area, and, for each sample point, shooting N_(T) raysat times t_(i,j) the time interval t₀ to t₀+T through the focus into thescene, with each ray spanning N_(L) sample points y_(i,j,k) on the lensarea A_(L). The manner in which subpixel jitter positions x_(i), pointsin time t_(i,j) and positions on the lens y_(i,j,k) are determined willbe described below. These three parameters determine the primary rayhitting the scene geometry in h(x_(i)t_(i,j),y_(i,j,k)) with the raydirection ω(x_(i),t_(i,j),y_(i,j,k)). In this manner, the value of theintegral in equation (7) can be approximated as

$\begin{matrix}{{L_{m,n} \approx {\frac{1}{N}{\sum\limits_{i = 0}^{N_{P} - 1}{\frac{1}{N_{T}}{\sum\limits_{j = 0}^{N_{T} - 1}{\frac{1}{N_{L}}{\sum\limits_{k = 0}^{N_{L} - 1}{{L\left( {{h\left( {x_{i},t_{i,j},y_{i,j,k}} \right)},{- {\omega \left( {x_{i},t_{i,j},y_{i,j,k}} \right)}}} \right)}{f_{m,n}\left( {x_{i},t_{i,j},y_{i,j,k}} \right)}}}}}}}}},} & (8)\end{matrix}$

where “N” is the total number of rays directed at the pixel.

It will be appreciated that rays directed from the scene toward theimage plane can comprise rays directly from one or more light sources inthe scene, as well as rays reflected off surfaces of objects in thescene. In addition, it will be appreciated that a ray that is reflectedoff a surface may have been directed to the surface directly from alight source, or a ray that was reflected off another surface. For asurface that reflects light rays, a reflection operator T_(fr) isdefined that includes a diffuse portion T_(fd), a glossy portion T_(fg)and a specular portion T_(fs), or

T _(ƒ) _(r) =T _(ƒd) +T _(ƒg) +T _(θs)  (9).

In that case, the Fredholm integral equation L=L_(e)+T_(fr)L governinglight transport can be represented as

L=L _(e) +T _(ƒ) _(r) _(−ƒ) _(s) L _(e) +T _(ƒ) _(s) (L−L _(e))+T _(ƒ)_(s) L+T _(ƒ) _(d) T _(ƒ) _(d) T _(ƒ) _(g+ƒ) _(s) L+T _(ƒ) _(d) T _(ƒ)_(d) L  (10),

where transparency has been ignored for the sake of simplicity;transparency is treated in an analogous manner. The individual terms inequation (10) are

(i) L_(e) represents flux due to a light source;

(ii) T_(ƒ) _(r) _(−ƒ) _(s) L_(e) (where T_(ƒ) _(r) _(−ƒ) _(s) =T_(ƒ)_(r) −T_(ƒ) _(s) ) represents direct illumination, that is, fluxreflected off a surface that was provided thereto directly by a lightsource; the specular component, associated with the specular portionT_(fs) of the reflection operator, will be treated separately since itis modeled using a δ-distribution;

(iii) T_(ƒ) _(g) (L−L_(e))representsglossyillumination,whichishandledbyrecursivedistribution raytracing, where, in the recursion, the source illumination has alreadybeen accounted for by the direct illumination (item (ii) above);

(iv) T_(ƒ) _(s) L represents a specular component, which is handled byrecursively using “L” for the reflected ray;

(v) T_(ƒ) _(d) T_(ƒ) _(g) _(+ƒ) _(s) L (where T_(ƒ) _(g) _(+ƒ) _(s)=T_(ƒ) _(g) +T_(ƒ) _(s) ) represents a caustic component, which is a raythat has been reflected off a glossy or specular surface (reference theT_(ƒ) _(g) _(+ƒ) _(s) operator) before hitting a diffuse surface(reference the T_(ƒ) _(d) operator); this contribution can beapproximated by a high resolution caustic photon map; and

(vi) T_(ƒ) _(d) T_(ƒ) _(d) L represents ambient light, which is verysmooth and is therefore approximated using a low resolution globalphoton map.

As noted above, the value of the integral (equation (7)) is approximatedby solving equation (8) making use of sample points x_(i), t_(i,j) andy_(i,j,k), where “x_(i)” refers to sample points within area A_(L) ofthe respective pixel at location (m,n) in the image plane, “t_(i,j)”refers to sample points within the time interval t₀ to t₀+T during whichthe shutter is open, and “y_(i,j,k)” refers to sample points on the lensA_(L). In accordance with one aspect of the invention, the sample pointsx_(i) comprise two-dimensional Hammersley points, which are defined as

$\left( {\frac{i}{N},{\Phi_{2}(i)}} \right),$

where 0≦i<N, and Φ₂(i) refers to the radical inverse of “i” inbase“two.” Generally, the “s” dimensional Hammersley point set is defined as

$\begin{matrix}{{\left. {{U_{N,s}^{Hammersley}\text{:}\mspace{14mu} \left\{ {0,\ldots \mspace{11mu},{N - 1}} \right\}}->I^{s}}i\mapsto x_{i} \right.:=\left( {\frac{i}{N},{\Phi_{b_{1}}(i)},\ldots \mspace{11mu},{\Phi_{b_{s - 1}}(i)}} \right)},} & (11)\end{matrix}$

where I^(s) is the s-dimensional unit cube [0,1)^(s) (that is, ans-dimensional cube each of whose dimension includes “zero,” and excludes“one”), the number of points “N” in the set is fixed and b₁, . . . ,b_(s−1), are bases. The bases do not need to be prime numbers, but theyare preferably relatively prime to provide a uniform distribution. Theradical inverse function Φ_(b), in turn, is generally defined as

$\begin{matrix}{{{\Phi_{b}\text{:}\mspace{14mu} N_{0}}->I}{i = \left. {\sum\limits_{j = 0}^{\infty}{{a_{j}(i)}b^{j}}}\mapsto{\sum\limits_{j = 0}^{\infty}{{a_{j}(i)}b^{{- j} - 1}}} \right.}} & (12)\end{matrix}$

where (a_(j))_(j=0) ^(∞) is the representation of “i” in integer base“b.” At N=(2^(n))², the two-dimensional Hammersley points are a (0, 2n,2)-net in base “two,” which are stratified on a 2^(n) by 2^(n) grid anda Latin hypercube sample at the same time. Considering the grid assubpixels, the complete subpixel grid underlying the image plane can befilled by simply abutting copies of the grid to each other.

Given integer subpixel coordinates (s_(x),s_(y)) the instance “i” andcoordinates (x,y) for the sample point x_(i) in the image plane can bedetermined as follows. Preliminarily, examining

$\left( {\frac{i}{N},{\Phi_{2}(i)}} \right)_{i = 0}^{N - 1},$

one observes that

(a) each line in the stratified pattern is a shifted copy of another,and

(b) the pattern is symmetric to the line y=x, that is, each column is ashifted copy of another column.

Accordingly, given the integer permutation σ(k):=2^(n)Φ₂(k) for0≦k<2^(n), subpixel coordinates (s_(x),s_(y)) are mapped onto stratacoordinates (j,k): (s_(x) mod 2^(n), s_(y) mod 2^(n)),an instance number“i” is computed as

i=j2^(n)+σ(k)  (13)

and the positions of the jittered subpixel sample points are determinedaccording to

$\begin{matrix}\begin{matrix}{x_{i} = \left( {{s_{x} + {\Phi_{2}(k)}},{s_{y} + {\Phi_{2}(j)}}} \right)} \\{= {\left( {{s_{x} + \frac{\sigma (k)}{2^{n}}},{s_{y} + \frac{\sigma (j)}{2^{n}}}} \right).}}\end{matrix} & (14)\end{matrix}$

An efficient algorithm for generating the positions of the jitteredsubpixel sample points x_(i) will be provided below in connection withCode Segment 1. A pattern of sample points whose positions aredetermined as described above in connection with equations (13) and (14)has an advantage of having much reduced discrepancy over a patterndetermined using a Halton sequence or windowed Halton sequence, asdescribed in the aforementioned Grabenstein application, and thereforethe approximation described above in connection with equation (8) givesin general a better estimation to the value of the integral describedabove in connection with equation (7). In addition, if “N” issufficiently large, sample points in adjacent pixels will have differentpatterns, reducing the likelihood that undesirable artifacts will begenerated in the image.

A ray tree is a collection of paths of light rays that are traced from apoint on the simulated camera's image plane into the scene. The computergraphics system 10 generates a ray tree by recursively followingtransmission, subsequent reflection and shadow rays using trajectorysplitting. In accordance with another aspect of the invention, a path isdetermined by the components of one vector of a global generalizedscrambled Hammersley point set. Generally, a scrambled Hammersley pointset reduces or eliminates a problem that can arise in connection withhigher-dimensioned low-discrepancy sequences since the radical inversefunction Φ_(b) typically has subsequences of b−1 equidistant valuesspaced by 1/b. Although these correlation patterns are merely noticeablein the full s-dimensional space, they are undesirable since they areprone to aliasing. The computer graphics system 10 attenuates thiseffect by scrambling, which corresponds to application of a permutationto the digits of the b-ary representation used in the radical inversion.For the symmetric permutation a from the symmetric group S_(b) overintegers 0, . . . , b−1, the scrambled radical inverse is defined as

$\begin{matrix}{{{\Phi_{b}\;:{N_{0} \times S_{b}}}->\left. {I\left( {i,\sigma} \right)}\mapsto{\sum\limits_{j = 0}^{\infty}{\sigma \; \left( {a_{j}(i)} \right)b^{{- j} - 1}}}\Leftrightarrow i \right.} = {\sum\limits_{j = 0}^{\infty}{{a_{j}(i)}{b^{j}.}}}} & (15)\end{matrix}$

If the symmetric permutation “σ” is the identity, the scrambled radicalinverse corresponds to the unscrambled radical inverse. In oneembodiment, the computer graphics system generates the symmetricpermutation a recursively as follows. Starting from the permutationσ₂=(0,1) for base b=2, the sequence of permutations is defined asfollows:

-   -   (i) if the base “b” is even, the permutation σ_(b) is generated        by first taking the values of

$2\sigma_{\frac{b}{2}}$

and appending the values of

${{2\sigma_{\frac{b}{2}}} + 1},$

an

-   -   (ii) if the base “b” is odd, the permutation σ_(b) is generated        by taking the values of σ_(b−1), incrementing each value that is        greater than or equal to

$\frac{b - 1}{2}$

by one, and inserting the value

$\frac{b - 1}{2}$

in the middle.This recursive procedure results in

σ₂=(0,1)

σ₃=(0,1,2)

σ₄=(0,2,1,3)

σ₅=(0,3,2,1,4)

σ₆=(0,2,4,1,3,5)

σ₇=(0,2,5,3,1,4,6)

σ₈=(0,4,2,6,1,5,3,7) . . .

The computer graphics system 10 can generate a generalizedlow-discrepancy point set as follows. It is often possible to obtain alow-discrepancy sequence by taking any rational s-dimensional point “x”as a starting point and determine a successor by applying thecorresponding incremental radical inverse function to the components of“x.” The result is referred to as the generalized low-discrepancy pointset. This can be applied to both the Halton sequence and the Hammersleysequence. In the case of the generalized Halton sequence, this can beformalized as

x _(i)=(Φ_(b) ₁ (i+i ₁), Φb ₂(i+i ₂), . . . , Φ(i+i _(s)))  (16),

where the integer vector (i₁, i₂, . . . , i_(s)) represents the offsetsper component and is fixed in advance for a generalized sequence. Theinteger vector can be determined by applying the inverse of the radicalinversion to the starting point “x.” A generalized Hammersley sequencecan be generated in an analogous manner.

Returning to trajectory splitting, generally trajectory splitting is theevaluation of a local integral, which is of small dimension and whichmakes the actual integrand smoother, which improves overall convergence.Applying replication, positions of low-discrepancy sample points aredetermined that can be used in evaluating the local integral. Thelow-discrepancy sample points are shifted by the corresponding elementsof the global scrambled Hammersley point set. Since trajectory splittingcan occur multiple times on the same level in the ray tree, branches ofthe ray tree are decorrelated in order to avoid artifacts, thedecorrelation being accomplished by generalizing the global scrambledHammersley point set.

An efficient algorithm for generating a ray tree will be provided belowin connection with Code Segment 2. Generally, in that algorithm, theinstance number “i” of the low-discrepancy vector, as determined abovein connection with equation (13), and the number “d” of used components,which corresponds to the current integral dimension, are added to thedata structure that is maintained for the respective ray in the raytree. The ray tree of a subpixel sample is completely specified by theinstance number “i.” After the dimension has been set to “two,” whichdetermines the component of the global Hammersley point set that is tobe used next, the primary ray is cast into the scene to span its raytree. In determining the determnuistic splitting by the components oflow discrepancy sample points, the computer graphics system 10 initiallyallocates the required number of dimensions “Δd.” For example, insimulating glossy scattering, the required number of dimensions willcorrespond to “two.” Thereafter, the computer graphics system 10generates scattering directions from the offset given by the scrambledradical inverses Φ_(b) _(d) (i,σ_(b) _(d) ), . . . , Φ_(b) _(d+Δd−1)(i,σ_(b) _(d+Δd−1) ), yielding the instances

$\begin{matrix}{{\left( y_{i,j} \right)_{j = 0}^{M - 1} = \left( {{{\Phi_{b_{d}}\left( {i,\sigma_{b_{d}}} \right)} \oplus \frac{j}{M}},\ldots \mspace{11mu},{{\Phi_{b_{d + {\Delta \; d} - 1}}\left( {i,\sigma_{b_{d + {\Delta \; d} - 1}}} \right)} \oplus {\Phi_{b_{d + {\Delta \; d} - 2}}\left( {j,\sigma_{b_{d + {\Delta \; d} - 2}}} \right)}}} \right)},} & (17)\end{matrix}$

where “⊕” refers to addition modulo “one.” Each direction of the “M”replicated rays is determined by y_(i,j) and enters the next level ofthe ray tree with d′:=d+Δd as the new integral dimension in order to usethe next elements of the low-discrepancy vector, and i′=i+j in order todecorrelate subsequent trajectories. Using an infinite sequenceoflow-discrepancy sample points, the replication heuristic is turnedinto an adaptive consistent sampling arrangement. That is, computergraphics system 10 can fix the sampling rate ΔM, compare current andprevious estimates every ΔM samples, and, if the estimates differ byless than a predetermined threshold value, terminate sampling. Thecomputer graphics system 10 can, in turn, determine the threshold value,by importance information, that is, how much the local integralcontributes to the global integral.

As noted above, the integral described above in connection with equation(7) is over a finite time period T from t₀ to t₀+T, during which timethe shutter of the simulated camera is open. During the time period, ifan object in the scene moves, the moving object may preferably bedepicted in the image as blurred, with the extent of blurring being afunction of the object's motion and the time interval t₀+T. Generally,motion during the time an image is recorded is linearly approximated bymotion vectors, in which case the integrand (equation (7)) is relativelysmooth over the time the shutter is open and is suited for correlatedsampling. For a ray instance “i,” started at the subpixel positionx_(i), the offset Φ₃(i) into the time interval is generated and theN_(T)−1 subsequent samples

${\Phi_{3}(i)} + \frac{j}{N_{T}}$

mod 1 are generated for 0<j<N_(T) that is

$\begin{matrix}{t_{i,j}:={t_{0} + {\left( {{\Phi_{3}(i)} \oplus \frac{j}{N_{T}}} \right) \cdot {T.}}}} & (18)\end{matrix}$

It will be appreciated that the value of N_(T) may be chosen to be“one,” in which case there will be no subsequent samples for rayinstance “i.” Determining sample points in this manner fills thesampling space, resulting in a more rapid convergence to the value ofthe integral (equation (7)). For subsequent trajectory splitting, raysare decorrelated by setting the instance i′=i+j.

In addition to determining the position of the jittered subpixel samplepoint x_(i), and adjusting the camera and scene according to the samplepoint t_(i,j) for the time, the computer graphics system also simulatesdepth of field. In simulating depth of field, the camera to be simulatedis assumed to be provided with a lens having known opticalcharacteristics and, using geometrical optics, the subpixel sample pointx_(i) is mapped through the lens to yield a mapped point x_(i)′. Thelens is sampled by mapping the dependent samples

$\begin{matrix}{y_{i,j,k} = \left( {\left( {{\Phi_{5}\left( {{i + j},\sigma_{5}} \right)} \oplus \frac{k}{N_{L}}} \right),\left( {{\Phi_{7}\left( {{i + j},\sigma_{7}} \right)} \oplus {\Phi_{2}(k)}} \right)} \right)} & (19)\end{matrix}$

onto the lens area A_(L) using a suitable one of a plurality of knowntransformations. As with N_(T), the value of N_(L) may be chosen to be“one.” Thereafter, a ray is shot from the sample point on the lensspecified by y_(i,j,k) through the point x_(i)′ into the scene. Theoffset (Φ₅(i+j,σ₅), Φ₇(i+j,σ₇)) in equation (19) comprises the nextcomponents taken from the generalized scrambled Hammersley point set,which, for trajectory splitting, is displaced by the elements

$\left( {\frac{k}{N_{L}},{\Phi_{2}(k)}} \right)$

of the two-dimensional Hammersley point set. The instance of the rayoriginating from sample point y_(i,j,k) is set to “i+j+k” in order todecorrelate further splitting down the ray tree. In equation (19), thescrambled samples (Φ₅(i+j,σ₅), Φ₇(1+j,σ₇)) are used instead of theunscrambled samples of (Φ₅(i+j), Φ₇(i+j)) since in bases “five” and“seven” up to five unscrambled samples will lie on a straight line,which will not be the case for the scrambled samples.

In connection with determination of a value for the direct illumination(T_(ƒ) _(r) _(−ƒ) _(s) L_(e) above), direct illumination is representedas an integral over the surface of the scene ∂V, which integral isdecomposed into a sum of integrals, each over one of the “L” single arealight sources in the scene. The individual integrals in turn areevaluated by dependent sampling, that is

$\begin{matrix}\begin{matrix}{{\left( {T_{{fr} - {fs}}L_{e}} \right)\left( {y,z} \right)} = {\int_{\partial V}{{L_{e}\left( {x,y} \right)}{f_{r}\left( {x,y,z} \right)}{G\left( {x,y} \right)}{x}}}} \\{= {\sum\limits_{k = 1}^{L}{\int_{{suppL}_{e,k}}{{L_{e}\left( {x,y} \right)}{f_{r}\left( {x,y,z} \right)}{G\left( {x,y} \right)}{x}}}}} \\{\approx {\sum\limits_{k = 1}^{L}{\frac{1}{M_{K}}{\sum\limits_{j = 0}^{M_{k} - 1}{{L_{e}\left( {x_{j},y} \right)}{f_{r}\left( {x_{j},y,z} \right)}{{G\left( {x_{j},y} \right)}.}}}}}}\end{matrix} & (20)\end{matrix}$

where suppL_(e,k) refers to the surface of the respective “k-th” lightsource. In evaluating the estimate of the integral for the “k-th” lightsource, for the M_(k)-th query ray, shadow rays determine the fractionof visibility of the area light source, since the point visibilityvaries much more than the smooth shadow effect. For each light source,the emission L_(e) is attenuated by a geometric term G, which includesthe visibility, and the surface properties are given by a bidirectionaldistribution function f_(r)−f_(s). These integrals are local integralsin the ray tree yielding the value of one node in the ray tree, and canbe efficiently evaluated using dependent sampling. In dependentsampling, the query ray comes with the next free integral dimension “d”and the instance “i,” from which the dependent samples are determined inaccordance with

$\begin{matrix}{x_{j} = {\left( {{{\Phi_{b_{d}}\left( {i,\sigma_{b_{d}}} \right)} \oplus \frac{j}{M_{k}}},{{\Phi_{b_{d + 1}}\left( {i,\sigma_{b_{d + 1}}} \right)} \oplus {\Phi_{2}(j)}}} \right).}} & (21)\end{matrix}$

The offset (Φ_(b) _(d) (i,σ_(b) _(d) ),Φ_(b) _(d+1) (i,σ_(b) _(d+1) ))again is taken from the corresponding generalized scrambled Hammersleypoint set, which shifts the two-dimensional Hammersley point set

$\left( {\frac{j}{M_{k}},{\Phi_{2}(j)}} \right)$

on the light source. Selecting the sample rate M_(k)=2^(n) ^(k) as apower of two, the local minima is obtained for the discrepancy of theHammersley point set that perfectly stratifies the light source. As analternative, the light source can be sampled using an arbitrarily-chosennumber M_(k) of sample points using

$\left( {\frac{j}{M_{k}},{\Phi_{M_{k}}\left( {j,\sigma_{M_{k}}} \right)}} \right)_{j - 0}^{M_{k} - 1}$

as a replication rule. Due to the implicit stratification of thepositions of the sample points as described above, the local convergencewill be very smooth.

The glossy contribution T_(ƒ) _(e) (L−L_(e)) is determined in a mannersimilar to that described above in connection with area light sources(equations (20) and (21)), except that a model f_(g) used to simulateglossy scattering will be used instead of the bidirectional distributionfunction f_(r). In determining the glossy contribution, two-dimensionalHammersley points are generated for a fixed splitting rate M and shiftedmodulo “one” by the offset (Φ_(b) _(d) (i,σ_(b) _(d) ), Φ_(b) _(d+1) ,(i,σ_(b) _(d+1) )) taken from the current ray tree depth given by thedimension field “d” of the incoming ray. The ray trees spanned into thescattering direction are decorrelated by assigning the instance fieldsi′=i+j in a manner similar to that done for simulation of motion blurand depth of field, as described above. The estimates generated for allrays are averaged by the splitting rate “M” and propagated up the raytree.

Volumetric effects are typically provided by performing a lineintegration along respective rays from their origins to the nearestsurface point in the scene. In providing for a volumetric effect, thecomputer graphics system 10 generates from the ray data a correspondingoffset Φ_(b) _(d) (i) which it then uses to shift the M equidistantsamples on the unit interval seen as a one-dimensional torus. In doingso, the rendering time is reduced in comparison to use of anuncorrelated jittering methodology. In addition, such equidistantshifted points typically obtain the best possible discrepancy in onedimension.

Global illumination includes a class of optical effects, such asindirect illumination, diffuse and glossy inter-reflections, causticsand color bleeding, that the computer graphics system 10 simulates ingenerating an image of objects in a scene. Simulation of globalillumination typically involves the evaluation of a rendering equation.For the general form of an illustrative rendering equation useful inglobal illumination simulation, namely:

L({right arrow over (x)},{right arrow over (w)})=L _(e)({right arrowover (x)},{right arrow over (w)})+∫_(s′)ƒ({right arrow over (x)},{rightarrow over (w)}′→{right arrow over (w)})G({right arrow over (x)},{rightarrow over (x)}′)V({right arrow over (x)},{right arrow over(x)}′)L({right arrow over (x)},{right arrow over (w)}′)dA′  (22)

it is recognized that the light radiated at a particular point {rightarrow over (x)} in a scene is generally the sum of two components,namely, the amount of light (if any) that is emitted from the point andthe amount of light (if any) that originates from all other points andwhich is reflected or otherwise scattered from the point {right arrowover (x)}. In equation (22), L({right arrow over (x)}, {right arrow over(w)}) represents the radiance at the point {right arrow over (x)} in thedirection {right arrow over (w)}=(θ,φ) (where “θ” represents the angleof direction {right arrow over (w)} relative to a direction orthogonalof the surface of the object in the scene containing the point {rightarrow over (x)}, and “φ” represents the angle of the component ofdirection {right arrow over (w)} in a plane tangential to the point{right arrow over (x)}). Similarly, L({right arrow over (x)}′,{rightarrow over (w)}′) in the integral represents the radiance at the point{right arrow over (x)}′ in the direction {right arrow over (w)}′=(θ′,φ′)(where “θ′” represents the angle of direction {right arrow over (w)}′relative to a direction orthogonal of the surface of the object in thescene containing the point {right arrow over (x)}′, and “φ′” representsthe angle of the component of direction {right arrow over (w)}′ in aplane tangential to the point {right arrow over (x)}′), and representsthe light, if any, that is emitted from point {right arrow over (x)}′which may be reflected or otherwise scattered from point {right arrowover (x)}.

Continuing with equation (22), L_(e)({right arrow over (x)},{right arrowover (w)}) represents the first component of the sum, namely, theradiance due to emission from the point {right arrow over (x)} in thedirection {right arrow over (w)}, and the integral over the sphere S′represents the second component, namely, the radiance due to scatteringof light at point {right arrow over (x)}. ƒ({right arrow over(x)},{right arrow over (w)}′→{right arrow over (w)}) is a bidirectionalscattering distribution function which describes how much of the lightcoming from direction {right arrow over (w)}′ is reflected, refracted orotherwise scattered in the direction {right arrow over (w)}, and isgenerally the sum of a diffuse component, a glossy component and aspecular component. In equation (22), the function G({right arrow over(x)}, {right arrow over (x)}′) is a geometric term

$\begin{matrix}{{G\left( {\overset{\rightarrow}{x},{\overset{\rightarrow}{x}}^{\prime}} \right)} = \frac{\cos \; {\theta cos\theta}^{\prime}}{{{\overset{\rightarrow}{x} - {\overset{\rightarrow}{x}}^{\prime}}}^{2}}} & (23)\end{matrix}$

where θ and θ′ are angles relative to the normals of the respectivesurfaces at points {right arrow over (x)} and {right arrow over (x)}′,respectively. Further in equation (22), V({right arrow over (x)}, {rightarrow over (x)}′) is a visibility function which equals the value one ifthe point {right arrow over (x)}′ is visible from the point {right arrowover (x)} and zero if the point {right arrow over (x)}′ is not visiblefrom the point {right arrow over (x)}.

The computer graphics system 10 makes use of global illuminationspecifically in connection with determination of the diffuse componentT_(ƒ) _(d) T_(ƒ) _(d) L and the caustic component T_(ƒ) _(d) T_(ƒ) _(g)_(+ƒ) _(s) L using a photon map technique. Generally, a photon map isconstructed by simulating the emission of photons by light source(s) inthe scene and tracing the path of each of the photons. For eachsimulated photon that strikes a surface of an object in the scene,information concerning the simulated photon is stored in a datastructure referred to as a photon map, including, for example, thesimulated photon's color, position and direction angle. Thereafter aRussian roulette operation is performed to determine the photon's state,that is, whether the simulated photon will be deemed to have beenabsorbed or reflected by the surface. If a simulated photon is deemed tohave been reflected by the surface, the simulated photon's direction isdetermined using, for example, a bidirectional reflectance distributionfunction (“BRDF”). If the reflected simulated photon strikes anothersurface, these operations will be repeated (reference the aforementionedGrabenstein application). The data structure in which information forthe various simulated photons is stored may have any convenient form;typically k-dimensional trees, for “k” an integer, are used. After thephoton map has been generated, it can be used in rendering therespective components of the image.

In generating a photon map, the computer graphics system 10 simulatesphoton trajectories, thus avoiding the necessity of discretizing thekernel of the underlying integral equation. The interactions of thephotons with the scene, as described above, are stored and used fordensity estimation. The computer graphics system 10 makes use of ascrambled low-discrepancy strictly-deterministic sequence, such as ascrambled Halton sequence, which has better discrepancy properties inhigher dimensions than does an unscrambled sequence. The scrambledsequence also has the benefit, over a random sequence, that theapproximation error decreases more smoothly, which will allow for use ofan adaptive termination scheme during generation of the estimate of theintegral. In addition, since the scrambled sequence is strictlydeterministic, generation of estimates can readily be parallelized byassigning certain segments of the low-discrepancy sequence to ones of aplurality of processors, which can operate on portions of thecomputation independently and in parallel. Since usually the space inwhich photons will be shot by selecting directions will be much largerthan the area of the light sources from which the photons were initiallyshot, it is advantageous to make use of components of smallerdiscrepancy, for example, Φ₂ or Φ₃(where, as above, Φ_(b) refers to theradical inverse function for base “b”), for use in connection withangles at which photons are shot, and components of higher discrepancy,for example, scrambled Φ₅ or Φ₇, for use in connection with sampling ofthe area of the respective light source, which will facilitate fillingthe space more uniformly.

The computer graphics system 10 estimates the radiance from the photonsin accordance with

$\begin{matrix}{{{\overset{\_}{L}}_{r}\left( {x,\omega} \right)} \approx {\frac{1}{A}{\sum\limits_{i \in {B_{k}{(x)}}}{{f_{r}\left( {\omega_{i},x,\omega} \right)}\Phi_{i}}}}} & (24)\end{matrix}$

where, in equation (24), Φ_(i) represents the energy of the respective“i-th” photon, ω_(i) is the direction of incidence of the “i-th photon,“B_(k)(x)” represents the set of the “k” nearest photons around thepoint “x,” and “A” represents an area around point “x” that includes thephotons in the set B_(k)(x). Since the energy of a photon is a functionof its wavelength, the Φ_(i) in equation (24) also represents the colorof the respective “i-th” photon. The computer graphics system 10 makesuse of an unbiased but consistent estimator for the area “A” for use inequation (24), which is determined as follows. Given a query ball, thatis, a sphere that is centered at point “x” and whose radius r(B_(k)(x))is the minimal radius necessary for the sphere to include the entire setB_(k)(x), a tangential disk D of radius r(B_(k)(x)) centered on thepoint “x” is divided into M equal-sized subdomains D_(i), that is

$\begin{matrix}{{{\bigcup_{i = 0}^{M - 1}D_{i}} = {{{D\mspace{14mu} {and}\mspace{14mu} D_{i}}\bigcap D_{j}} \neq {0\mspace{14mu} {for}\mspace{14mu} i} \neq j}},\mspace{14mu} {{{where}\mspace{14mu} {D_{i}}} = {\frac{D}{M} = {\frac{\pi \; {r^{2}\left( {B_{k}(x)} \right)}}{M}.}}}} & (25)\end{matrix}$

The set

P={D _(i) |D _(i) ∩{x _(i)|_(D) |iεB _(k)(x)}≠0}  (26)

contains all the subdomains D_(i) that contain a point x_(i)|_(D) on thedisk, which is the position of the “i-th” photon projected onto theplane defined by the disk D along its direction of incidence ω_(i).Preferably, the number M of subdomains will be on the order of √{squareroot over (k)} and the angular subdivision will be finer than the radialsubdivision in order to capture geometry borders. The actual area A isthen determined by

$\begin{matrix}{A = {\pi \; {r^{2}\left( {B_{k}(x)} \right)}{\frac{P}{M}.}}} & (27)\end{matrix}$

Determining the actual coverage of the disk D by photons significantlyimproves the radiance estimate (equation (24)) in corners and onborders, where the area obtained by the standard estimate πr²(B_(k)(x))would be too small, which would be the case at corners, or too large,which would be the case at borders. In order to avoid blurring of sharpcontours of caustics and shadows, the computer graphics system 10 setsthe radiance estimate L to black if all domains D_(i) that touch point“x” do not contain any photons.

It will be appreciated that, in regions of high photon density, the “k”nearest photons may lead to a radius r(B_(k)(x) that is nearly zero,which can cause an over-modulation of the estimate. Over-modulation canbe avoided by selecting a minimum radius r_(min), which will be used ifr(B_(k)(x)) is less than r_(min). In that case, instead of equation(24), the estimate is generated in accordance with

$\begin{matrix}{{{\overset{\_}{L}}_{r}\left( {x,\omega} \right)} = {\frac{N}{k}{\sum\limits_{i \in {B_{k}{(x)}}}{\Phi_{i}{f_{r}\left( {\omega_{i},x,\omega_{r}} \right)}}}}} & (28)\end{matrix}$

assuming each photon is started with 1/N of the total flux Φ. Theestimator in equation (28) provides an estimate for the mean flux of the“k” photons if r(B_(k)(x))<r_(min).

The global photon map is generally rather coarse and, as a result,subpixel samples can result in identical photon map queries. As aresult, the direct visualization of the global photon map is blurry andit is advantageous to perform a smoothing operation in connectiontherewith In performing such an operation, the computer graphics system10 performs a local pass integration that removes artifacts of thedirect visualization. Accordingly, the computer graphics system 10generates an approximation for the diffuse illumination term T_(ƒ) _(d)T_(ƒ) _(d) L as

$\begin{matrix}{{{{T_{f_{d}}T_{f_{d}}L} \approx {\left( {T_{f_{d}}{\overset{\_}{L}}_{r}} \right)(x)}} = {{\int_{S^{2}{(x)}}{{f_{d}\ (x)}{{\overset{\_}{L}}_{r}\left( {h\left( {x,\overset{\rightarrow}{\omega}} \right)} \right)}\cos \; \theta \; {\omega}}} \approx {\frac{f_{d}(x)}{M}{\sum\limits_{i = 0}^{M - 1}{{\overset{\_}{L}}_{r}\left( {h\left( {x,{\omega \left( {{\arcsin \sqrt{u_{i,1}}},{2\pi \; u_{i,2}}} \right)}} \right)} \right)}}}}},} & (29)\end{matrix}$

with the integral over the hemisphere S²(x) of incident directionsaligned by the surface normal in “x” being evaluated using importancesampling. The computer graphics system 10 stratifies the sample pointson a two-dimensional grid by applying dependent trajectory splittingwith the Hammersley sequence and thereafter applies irradianceinterpolation. Instead of storing the incident flux Φ_(i) of therespective photons, the computer graphics system 10 stores theirreflected diffuse power f_(d)(x_(i))Φ_(i) with the respective photons inthe photon map, which allows for a more exact approximation than can beobtained by only sampling the diffuse BRDF in the hit points of thefinal gather rays. In addition, the BRDF evaluation is needed only onceper photon, saving the evaluations during the final gathering. Insteadof sampling the full grid, the computer graphics system 10 uses adaptivesampling, in which refinement is triggered by contrast, distancetraveled by the final gather rays in order to more evenly sample theprojected solid angle, and the number of photons that are incident formthe portion of the projected hemisphere. The computer graphics systemfills in positions in the grid that are not sampled by interpolation.The resulting image matrix of the projected hemisphere is medianfiltered in order to remove weak singularities, after which theapproximation is generated. The computer graphics system 10 performs thesame operation in connection with, for example, hemispherical skyillumination, spherical high dynamic-range environmental maps, or anyother environmental light source.

The computer graphics system 10 processes final gather rays that strikeobjects that do not cause caustics, such as plane glass windows, byrecursive ray tracing. If the hit point of a final gather ray is closerto its origin than a predetermined threshold, the computer graphicssystem 10 also performs recursive ray tracing. This reduces thelikelihood that blurry artifacts will appear in corners, which mightotherwise occur since for close hit points the same photons would becollected, which can indirectly expose the blurry structure of theglobal photon map.

Generally, photon maps have been taken as a snapshot at one point intime, and thus were unsuitable in connection with rendering of motionblur. Following the observation that averaging the result of a pluralityof photon maps is generally similar to querying one photon map with thetotal number of photons from all of the plurality of photon maps, thecomputer graphics system 10 generates N_(T) photon maps, where N_(T) isdetermined as described above, at points in time

$\begin{matrix}{{t_{b} = {t_{0} + {\frac{b + \frac{1}{2}}{N_{T}}T}}},} & (30)\end{matrix}$

for 0≦b<N_(T). As noted above, N_(T) can equal “one,” in which case “N”photon maps are used, with “N” being chosen as described above. In thatcase,

t _(i) =t ₀+Φ₃(i)T  (31),

and thus t_(i,j)=t_(i,0), that is, “t_(i),” for N_(T)=1. In the generalcase (equation (30)), during rendering, the computer graphics system 10uses the photon map with the smallest time difference |t_(i,j)−t_(b)| inconnection with rendering for the time sample point t_(i,j).

The invention provides a number of advantages. In particular, theinvention provides a computer graphics system that makes use of strictlydeterministic distributed ray tracing based on low-discrepancy samplingand dependent trajectory splitting in connection with rendering of animage of a scene. Generally, strictly deterministic distributed raytracing based on deterministic low-discrepancy sampling and dependenttrajectory splitting is simpler to implement than an implementationbased on random or pseudo-random numbers. Due to the properties of theradical inverse function, stratification of sample points is intrinsicand does not need to be considered independently of the generation ofthe positions of the sample points. In addition, since the methodologyis strictly deterministic, it can be readily parallelized by dividingthe image into a plurality of tasks, which can be executed by aplurality of processors in parallel. There is no need to take a step ofensuring that positions of sample points are not correlated, which isgenerally necessary if a methodology based on random or pseudo-randomnumbers is to be implemented for processing in parallel.

Moreover, the methodology can be readily implemented in hardware, suchas a graphics accelerator, particularly if Hammersley point sets areused, since all points with a fixed index “i” yield a regular grid. Agraphics accelerator can render a plurality of partial imagescorresponding to these regular grids in a number of parallel tasks, andinterleave the partial images in an accumulation buffer to provide thefinal image. Operating in this manner provides very good load balancingamong the parallel tasks, since all of the tasks render almost the sameimage.

In addition, the methodology can readily be extended to facilitaterendering of animations. Generally, an animation consists of a series offrames, each frame comprising an image. In order to decorrelate thevarious frames, instead of initializing the field of integers used asidentifiers for ray instances for each frame by “i,” “i+i_(f)” can beused, where “i_(f)” is a frame number. This operates as an offsetting of“i” by “i_(f)” which is simply a generalization of the Hammersleypoints. A user can select to initialize the field of integers for eachframe by “i,” in which case the frames will not be correlated. In thatcase, undersampling artifacts caused by smooth motion will remain localand are only smoothly varying. Alternatively, the user can select toinitialize the field of integers for each frame by “i+i_(f)” in whichcase the artifacts will not remain local, and will instead appear asnoise or film grain flicker in the final animation. The latter issometimes a desired feature of the resulting animation, whether forartistic reasons or to match actual film grain. Another variation is toadd if directly to k and clip the result by 2^(n) (reference CodeSegment 1, below). In that case, the pixel sampling pattern will changefrom frame to frame and the frame number i_(f) will need to be known inthe post-production process in order to reconstruct the pixel samplingpattern for compositing purposes.

Generally, a computer graphics system that makes use of deterministiclow-discrepancy sampling in determination of sample points will performbetter than a computer graphics system that makes use of random orpseudo-random sampling, but the performance may degrade to that of asystem that makes use of random or pseudo-random sampling in higherdimensions. By providing that the computer graphics system performsdependent splitting by replication, the superior convergence oflow-dimensional low-discrepancy sampling can be exploited with theeffect that the overall integrand becomes smoother, resulting in betterconvergence than with stratified random or pseudo-random sampling. Sincethe computer graphics system also makes use of dependent trajectorysampling by means of infinite low discrepancy sequences, consistentadaptive sampling of, for example, light sources, can also be performed.

In addition, it will be appreciated that, although the computer graphicssystem has been described as making use of sample points generated usinggeneralized scrambled and/or unscrambled Hammersley and Haltonsequences, it will be appreciated that generally any (t,m,s)-net or(t,s)-sequence can be used.

At a more general level, the invention provides an improved quasi-MonteCarlo methodology for evaluating an integral of a function ƒ on the “s”dimensional unit cube [0,1)^(s). In contrast with this methodology,which will be referred to as trajectory splitting by dependentsplitting, in prior methodologies, the sample points in the integrationdomain for which the sample values at which sample values for thefunction were generated were determined by providing the same number ofcoordinate samples along each dimension. However, for some dimensions ofan integrand, it is often the case that the function ƒ will exhibit ahigher variance than for other dimensions. The invention exploits thisby making use of trajectory splitting by dependent samples in criticalregions.

For example, assume that a function ƒ is to be integrated over s=s₁+s₂dimensions. The partial integral (equation (32))

$\begin{matrix}{{g(x)} = {{\int_{I^{s_{2}}}{{f\left( {x,y} \right)}\ {y}}} \approx {\frac{1}{N_{2}}{\sum\limits_{j = 0}^{N_{2} - 1}{f\left( {x,y_{j}} \right)}}}}} & (32)\end{matrix}$

(x″ and “y” comprising disjoint sets of the “s” dimensions, and x∪ycomprising the set of all of the dimensions), where “N₂” identifies thenumber of samples selected for the set “y” of dimensions, can be definedover the portion of the integration domain that is defined by unit cube[0,i)^(s) ² , which, in turn, corresponds to the portion of theintegration domain that is associated with set s₂ dimensions. Evaluatingg(x) using equation (32) will affect a smoothing of the function ƒ inthe s₂ dimensions that are associated with set “y.”

The result generated by applying equation (32) can then be used toevaluate the full integral

$\begin{matrix}{{{\int_{I^{s_{1}}}{\int_{I^{s_{2}}}{{f\left( {x,y} \right)}\ {y}\mspace{11mu} {x}}}} = {{\int_{I^{s_{1}}}{{g(x)}\ {x}}} \approx {\frac{1}{N_{1}}{\sum\limits_{i = 0}^{N_{1} - 1}{\frac{1}{N_{2}}{\sum\limits_{j = 0}^{N_{2} - 1}{f\left( {x_{i},y_{j}} \right)}}}}}}},} & (33)\end{matrix}$

where “N₁” identifies the number of samples selected for the set “x” ofdimensions, that is, over the remaining dimensions of the integrationdomain that are associated with the s₁ dimensions that are associatedwith the set “x.” If the dimension splitting x, y is selected such thatthe function ƒ exhibits relatively high variance over the set “y” of theintegration domain, and relatively low variance over the set “x,” itwill not be necessary to generate sample values for the functionN₁-times-N₂ times. In that case, it will suffice to generate sample onlyvalues N₂ times over the integration domain. If the correlationcoefficient of ƒ(ξ,η) and ƒ(ξ,η′), which indicates the degree ofcorrelation between values of function evaluated, for the former, at(x_(i),y_(i))=(ξ,η), and, for the later, at (x_(i),y_(i))=(ξ,η′), isrelatively high, the time complexity required to evaluate the functionƒ_([0,1)) _(s) (x₀, . . . , x_(s−1)) will be decreased.

The smoothness of an integrand can be exploited using a methodology thatwill be referred to as correlated sampling. Generally, that is, ifcorrelated sampling is not used, in evaluating an integral, eachdimension will be associated with its respective sequence. However, incorrelated sampling, the same sequence can be used for all of thedimensions over the integration domain, that is

$\begin{matrix}{{{\frac{1}{M}{\sum\limits_{j = 1}^{N}{\frac{1}{N_{j}}{\sum\limits_{i = 0}^{N_{j} - 1}{f_{j}\left( {x_{i},y_{j}} \right)}}}}} \approx {\frac{1}{M}{\sum\limits_{j = 1}^{M}{\int_{I^{s}}{{f_{j}(x)}\ {x}}}}}} = {{\int_{I^{s}}{\frac{1}{M}{\sum\limits_{j = 1}^{M}{{f_{j}(x)}\ {x}}}}} \approx {\frac{1}{N}{\sum\limits_{i = 0}^{N - 1}{\frac{1}{M}{\sum\limits_{j = 1}^{M}{f_{j}\left( x_{i} \right)}}}}}}} & (34)\end{matrix}$

The methodology of trajectory splitting by depending sampling makes useof a combination of the trajectory splitting technique described abovein connection with equations (32) and (33) with the correlated samplingmethodology described in connection with equation (34).

Since integrals are invariant under toroidal shifting for z_(j)εI^(s) ², that is,

$\begin{matrix}{{\left. \begin{matrix}{S_{j}:\left. I^{s_{2}}\rightarrow I^{s_{2}} \right.} \\\left. y\mapsto{\left( {y + z_{j}} \right)\mspace{11mu} {mod}\; 1} \right.\end{matrix}\Rightarrow{\int_{I^{s_{2}}}{{g(y)}\ {y}}} \right. = {\int_{I^{s_{2}}}{{g\left( {S_{j}(y)} \right)}\ {y}}}},} & (35)\end{matrix}$

the values of the integrals also do not change. Thus, if, in equation(33), the inner integral is replicated “M” times,

$\begin{matrix}\begin{matrix}{{\int_{I^{s_{1}}}{\int_{I^{s_{2}}}{{f\left( {x,y} \right)}\ {y}\mspace{11mu} {x}}}} = {\int_{I^{s_{1}}}{\int_{I^{s_{2}}\;}{\frac{1}{M}\ {\sum\limits_{j = 0}^{M - 1}{{f\left( {x,{S_{j}(y)}} \right)}{y}\mspace{11mu} {x}}}}}}} \\{\approx {\frac{1}{N}{\sum\limits_{i = 0}^{N - 1}{\frac{1}{M}{\sum\limits_{j = 0}^{M - 1}{f\left( {x_{i},{S_{j}\left( y_{i} \right)}} \right)}}}}}} \\{= {\frac{1}{N}{\sum\limits_{i = 0}^{N - 1}{\frac{1}{M}{\sum\limits_{j = 0}^{M - 1}{{f\left( {x_{i},{\left( {y_{i} + z_{j}} \right)\mspace{11mu} {mod}\; 1}} \right)}.}}}}}}\end{matrix} & (36)\end{matrix}$

For index “j,” the functions ƒ(x_(i), S_(j)(y_(i))) are correlated,enabling the smoothness of the integrand in those dimensions that arerepresented by “y” to be exploited, as illustrated above in connectionwith equation (19) (lens sampling), equations (20) and (21) (area lightsources) and equation (29) (approximation for the diffuse illuminationterm). It will be appreciated that the evaluation using the replicationis the repeated application of the local quadrature rule U_(M,s) ₂:=(z_(j))_(j=0) ^(M) shifted by random offset values y_(i). The use ofdependent variables in this manner pays off particularly if there issome smoothness in the integrand along one or more dimensions. Splittingcan be applied recursively, which yields a history tree, in which eachpath through the respective history tree represents a trajectory of aparticle such as a photon.

The quasi-Monte Carlo methodology of trajectory splitting by dependentsampling makes use of sets of deterministic, low-discrepancy samplepoints both for the global quadrature rule U_(N,s) ₁ _(+s) ₂=(x_(i),y_(i))_(i=0) ^(N), that is, integration over all of thedimensions s₁+s₂ comprising the entire integration domain, as well asfor the local quadrature rule U_(M,s) ₂ :=(z_(j))_(j=0) ^(M), that is,integration over the dimensions s₂ of the integration domain. Themethodology unites splitting and dependent sampling, exploiting thestratification properties of low-discrepancy sampling. Accordingly, itwill be possible to concentrate more samples along those dimensions inwhich the integrand exhibits high levels of variation, and fewer samplesalong those dimensions in which the integrand exhibits low levels ofvariation, which reduces the number of sample points at which thefunction will need to be evaluated. If the methodology is to be appliedrecursively a plurality of times, it will generally be worthwhile tocalculate a series of values z_(j) that are to comprise the set U_(M,s)₂ . In addition, the methodology may be used along with importancesampling and, if U is an infinite sequence, adaptive sampling. Inconnection with adaptive sampling, the adaptations will be applied inthe replication independently of the sampling rate, so that thealgorithm will remain consistent. The low-discrepancy sample points setsU_(N,s) ₁ _(+s) ₂ and U_(M,s) ₂ can be chosen arbitrarily; for example,the sample point set U_(M,s) ₂ can be a projection of sample point setU_(N,s) ₁ _(+s) ₂ . When trajectory splitting is recursively applied tobuild trajectory trees, generalizing the point set U_(N,s) ₁ _(+s) ₂ forthe subsequent branches can be used to decorrelate the separate parts ofthe respective tree.

Code Segment 1

The following is a code fragment in the C++ programming language forgenerating the positions of the jittered subpixel sample points x_(i)

unsigned short Period, *Sigma; void InitSigma(int n) { unsigned shortInverse, Digit, Bits; Period = 1 << n; Sigma = new unsigned short[Period]; for (unsigned short i = 0; i < Period; i++) { Digit = PeriodInverse = 0; for (bits = i; bits; bits >>= 1) { Digit >>= 1; if (Bits&1) inverse += Digit; } Sigma[i] = Inverse;  } } voidSampleSubpixel(unsigned int *i, double *x, double *y, int s_(x), ints_(y)) { int j = s_(x) & (Period − 1); int k = s_(y) & (Period − 1); *i= j * Period + Sigma[k] *x = (double) s_(x) + (double) Sigma[k] /(double) Period; *y = (double) s_(y) + (double) Sigma[j] / (double)Period; }

Code Segment 2

The following is a code fragment in the C++ programming language forgenerating a ray tree

class Ray { int i; //current instance of low discrepancy vector int d;//current integral dimension in ray tree . . . } void Shade (Ray & Ray){ Ray next_ray; int i = ray.i; int d = ray.d . . . for (int j = 0; j <M; j++) { . . . //ray set up for recursion$y = \left( {{{\Phi_{b_{d}}\left( {i,\sigma_{b_{d}}} \right)} \oplus \frac{j}{M}},\ldots \mspace{14mu},{{\Phi_{b_{d + {\Delta \; d} - 1}}\left( {i,\sigma_{b_{d + {\Delta \; d} - 1}}} \right)} \oplus {\Phi_{b_{d + {\Delta \; d} - 2}}\left( {i,\sigma_{b_{d + {\Delta \; d} - 2}}} \right)}}} \right)$next_ray = SetUpRay(y); //determine ray parameters by y next_ray.i = i +j; //decorrelation by generalization next_ray.d = d + Δd; //dimensionallocation Shade(next_ray); . . . } }

It will be appreciated that a system in accordance with the inventioncan be constructed in whole or in part from special purpose hardware ora general purpose computer system, or any combination thereof, anyportion of which may be controlled by a suitable program. Any programmay in whole or in part comprise part of or be stored on the system in aconventional manner, or it may in whole or in part be provided in to thesystem over a network or other mechanism for transferring information ina conventional manner. In addition, it will be appreciated that thesystem may be operated and/or otherwise controlled by means ofinformation provided by an operator using operator input elements (notshown) which may be connected directly to the system or which maytransfer the information to the system over a network or other mechanismfor transferring information in a conventional manner.

The foregoing description has been limited to a specific embodiment ofthis invention. It will be apparent, however, that various variationsand modifications may be made to the invention, with the attainment ofsome or all of the advantages of the invention. It is the object of theappended claims to cover these and such other variations andmodifications as come within the true spirit and scope of the invention.

1. A computer graphics system for generating a pixel value for a pixelin an image, the pixel value being representative of a point in a sceneas recorded on an image plane of a simulated camera, the computergraphics system comprising: A. a sample point generator configured togenerate a set of sample points, at least one sample points beinggenerated using at least one dependent sample, the at least onedependent sample comprising at least one element of a low-discrepancysequence offset by at least one element of another low-discrepancysequence; and B. a function evaluator configured to generate at leastone value representing an evaluation of a selected function at one ofthe sample points generated by said sample point generator, the valuegenerated by the function evaluator corresponding to the pixel value. 2.A computer-implemented method of generating a pixel value for a pixel inan image, the pixel value being representative of a point in a scene asrecorded on an image plane of a simulated camera, the method comprising:A. a sample point generation step of generating a set of sample points,at least one sample points being generated using at least one dependentsample, the at least one dependent sample comprising at least oneelement of a low-discrepancy sequence offset by at least one element ofanother low-discrepancy sequence; and B. a function evaluation step ofgenerating at least one value representing an evaluation of a selectedfunction at one of the sample points generated during the sample pointgeneration step, the value generated during the function evaluation stepcorresponding to the pixel value.
 3. A computer program product for usein connection with a computer to provide a computer graphics system forgenerating a pixel value for a pixel in an image, the pixel value beingrepresentative of a point in a scene as recorded on an image plane of asimulated camera, the computer program product comprising acomputer-readable medium having encoded thereon: A. a sample pointgenerator module configured to enable the computer to generate a set ofsample points, at least one sample points being generated using at leastone dependent sample, the at least one dependent sample comprising atleast one element of a low-discrepancy sequence offset by at least oneelement of another low-discrepancy sequence; and B. a function evaluatormodule configured to enable the computer to generate at least one valuerepresenting an evaluation of a selected function at one of the samplepoints generated by said sample point generator, the value generated bythe computer corresponding to the pixel value.